Simple denoising algorithm using wavelet 

transform 



Manojit Roy, V. Ravi Kumar , B. D. Kulkarni 

Chemical Engineering Division, National Chemical Laboratory, 

Pune 411 008, India 
John Sanderson, Martin Rhodes 
Department of Chemical Engineering, Monash University, 
Clayton, Victoria, 3168, Australia 
Michel vander Stappen 
Unilever Research, Vlaardingen, Postbox 114, 3130, 
AC Vlaardingen, The Netherlands 

February 6, 2008 



keywords: Noise, Discrete wavelet transform, Chaos, Differentiation 

Application of wavelets and multiresolution analysis to reaction engineering 
systems from the point of view of process monitoring, fault detection, systems 
analysis etc. is an important topic and of current research interest (see, Bakshi 
and Stephanopoulos, 1994; Safavi et. a!., 1997; Luo et. al., 1998; Carrier and 
Stephanopoulos, 1998). In the present paper we focus on one such important 
application, where we propose a new and simple algorithm for the reduction 
of noise from a scalar time series data. Presence of noise in a time-varying 
signal restricts one's ability to obtain meaningful information from the sig- 
nal. Measurement of correlation dimension can get affected by a noise level as 
small as 1% of signal, making estimation of invariant properties of a dynamical 
^-mail address for correspondence: ravi@che.ncl. res. in 
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system, such as the dimension of the attractor and Lyapunov exponents, al- 
most impossible (Kostelich and Yorke, 1988). Noise in experimental data can 
also cause misleading conclusions (Grassberger et. al., 1991). A host of lit- 
erature exists on various techniques for noise reduction (Kostelich and Yorke, 
1988; Hardle, 1990; Farmer and Sidorowich, 1991; Sauer, 1992; Cawley and 
Hsu, 1992; Cohen, 1995; Donoho and Johnstone, 1995; Kantz and Schreiber, 
1997). For instance, Fast Fourier Transform (FFT) reduces noise effectively in 
those cases where the frequency distribution of noise is known (Kostelich and 
Yorke, 1988; Cohen, 1995; Kantz and Schreiber, 1997); singular value analy- 
sis methods (Cawley and Hsu, 1992) project the original time-series onto an 
optimal subspace, whereby noise components are left behind in the remaining 
orthogonal directions, etc. In the existing wavelet-based denoising methods 
(Donoho and Johnstone, 1995) two types of denoising are introduced: linear 
denoising and nonlinear denoising. In linear denoising, noise is assumed to 
be concentrated only on the fine scales and all the wavelet coefficients below 
these scales are cut off. Nonlinear denoising, on the other hand, treats noise 
reduction by either cutting off all coefficients below a certain threshold (so 
called 'hard-thresholding'), or reducing all coefficients by this threshold (so 
called 'soft-thresholding'). The threshold values are obtained by statistical 
calculations and has been seen to depend on the standard deviation of the 
noise (Nason, 1994). 

The noise reduction algorithm that we propose here makes use of the 
wavelet transform (WT) which in many ways complements the well known 
Fourier Transform (FT) procedure. We apply our method, firstly, to three 
model flow systems, viz. Lorenz, Autocatalator, and Rossler systems, all ex- 
hibiting chaotic dynamics. The reasons for choosing these systems are the 
following: Firstly, all of them are simplified models of well-studied experimen- 
tal systems. For instance, Lorenz is a simple realization of convective systems 
(Lorenz, 1963), while the Autocatalator and Rossler have their more compli- 
cated analogs in chemical multicomponent reactions (Rossler, 1976; Lynch, 
1992). Secondly, chaotic dynamics is extremely nonlinear, highly sensitive, 
possesses only short-time correlations and is associated with a broad range 
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of frequencies (Guckenheimer and Holmes, 1983; Strogatz, 1994). Because of 
these properties it is well known that FT methods are not applicable in a 
straightforward way to chaotic dynamical systems (Abarbanel, 1993). On the 
other hand, WT methods are particularly suited to handle not only nonlin- 
ear but nonstationary signals (Strang and Nguyen, 1996). This is because the 
properties of the data are studied at varying scales with superior time localiza- 
tion analysis when compared to FT technique. Our noise reduction algorithm 
is advantageous, because, as shall be shown, the threshold level for noise is 
identified automatically. In this study, we have used the discrete analog of 
the wavelet transform (DWT) which involves transforming a given signal with 
orthogonal wavelet basis functions by dilating and translating in discrete steps 
(Daubechies, 1990; Holschneider, 1995). For study purposes we corrupt one 
variable x(t) for each of these systems with noise of zero mean, and then apply 
our algorithm for denoising. We analyze the performance of this method in all 
the three systems for a wide range of noise strengths, and show its effectiveness. 
Importantly, we then validate the applicability of the method to experimental 
data obtained from two chemical systems. In one system the time series data 
was obtained from pressure fluctuation measurements of the hydrodynamics 
in a fluidized bed. In the other the conductivity measurements in a liquid 
surfactant manufacturing experiment were analyzed. 

Methodology 

The noise reduction algorithm based on DWT consists of the following five 
steps: 

Step 1: 

In first step, we differentiate the noisy signal x(t) to obtain the data Xd(i), 
using the method of central finite differences with fourth order correction to 
minimize error (Constantinides, 1987), i.e., 

*(«> = m 
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Step 2: 



We then take DWT of the data xa(t) and obtain wavelet coefficients W^ k 
at various dyadic scales j and displacements k. A dyadic scale is the scale 
whose numerical magnitude is equal to 2 (two) raised to an integer exponent, 
and is labeled by the exponent. Thus, the dyadic scale j refers to a scale of 
magnitude 2 j . In other words, it indicates a resolution of 2 j data points. Thus 
a low value of j implies finer resolution while high j analyzes the signal at 
coarser resolution. This transform is the discrete analog of continuous WT 
(Holschneider, 1995), and is given by the formula 

/+oo 
x d (t) if> jtk (t) dt , (2) 
-oo 

with 

^■ fc (f) = 2^(2H-k) 

where j, k are integers. As for the wavelet function ip{t) we have chosen 
Daubechies compactly supported orthogonal function with four filter coeffi- 
cients (Daubechies, 1990; Press et. al., 1996). 

Step 3: 

In this step we estimate the power Pj contained in different dyadic scales j, 
via 

+oo 

Pj {x) = £ \W jtk \ 2 (j = l,2,-..) (3) 

k=— oo 

By plotting the variation of Pj with j, we see that it is possible to identify a 
scale j m at which the power due to noise falls off rapidly. This is important 
because as we shall see from the studies of the case examples that it provides a 
mean for automation in detection of threshold. The identification of the scale 
j m at which power due to noise shows the first minimum allows us to reset all 
W jtk upto scale j m to zero, i.e., W jtk = 0, for j = 1, 2, • • -,j m . 
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Step 4 : 



In the fourth step, we reconstruct the denoised data x d {t) by taking inverse 
transform of the coefficients Wj t t : 

oo oo 

Xd(t) = E w i,k *Pj,k(t) , (4) 

j=0 fe=— oo 

where is normalization constant given by 

H = VS-J^du < oo, 

with ^(u) as the Fourier transform of the wavelet function ip(t). 

Step 5: 

In the fifth and final step x d (t) is integrated to yield the cleansed signal x(t): 

x(t) = J x d (t)dt . (5) 

There exists a commutativity property between the operation of differen- 
tiation/integration and wavelet transform. Therefore first differentiating the 
signal and then taking DWT is equivalent to carrying out the two operations in 
reverse order. This implies that the same result can be obtained by switching 
the order between the first and second steps, and then between the fourth and 
fifth. 

The effectiveness of the method lies in the following observations. Upon 
differentiation, contribution due to white noise moves towards the finer scales 
because the process of differentiation converts the uncorrelated stochastic pro- 
cess to a first order moving average process and thereby distributes more energy 
to the finer scales. That the differentiation of white noise brings about this be- 
havior is known in the Fourier spectrum (Box et. al., 1994). It may be noted 
that the nature and effectiveness of separation depend on the wavelet basis 
function chosen and also on the properties of the derivatives of WT, which 
is in itself a highly interesting and not fully understood subject (Strang and 
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Nguyen, 1996). For the signals studied in this paper, model as well as exper- 
imental, the wanted signal features lie in the coarser wavelet scales while the 
unwanted signal features after differentiation lie in the finer resolution wavelet 
scales. This is because the size of the data set handled decides the total num- 
ber of scales available and a suitable choice can bring out the noisy signal WT 
features lying in the coarser scales. This justifies the assumption that fine scale 
features can be removed by setting the corresponding wavelet coefficients to 
zero and coarse scale features retained after differentiation. For this reason we 
also see a clear separation in the scales attributed to noise and those for the 
signal. A threshold scale for noise removal is thus identified and this leads to 
an automation for noise removal. 

Result and Discussion 

We first take up the three model systems and discuss the observations. For 
test purposes the pure signal obtained from these systems were corrupted with 
noise of certain strength. For the systems chosen for study, viz., Lorenz, Au- 
tocatalator, and Rossler, Table [I] summarizes the details, i.e. the equations 
governing their dynamics, the values chosen for the parameters, and the na- 
ture of the evolution of these systems for these sets of parameter values. These 
values were chosen appropriately so that the dynamics is chaotic. In our ini- 
tial studies, purely for testing purposes, we studied situations where we ensure 
that all scales are affected by noise. In the wavelet domain this can be conve- 
niently carried out by perturbing the wavelet coefficients in the following way. 
The differential equations are first numerically integrated to obtain pure signal 
x°(U) at equidistant time steps U. We then take DWT of the signal, and add 
white noise rj of zero mean and certain strength, i.e., Wj t k = Wj k + rj, where 
Wj k , Wj t k are the wavelet coefficients of pure and noisy signals respectively. 
We take the strength of noise as the relative percentage of the difference be- 
tween the maximum and minimum of the signal value. Since each coefficient 
Wj t k is individually affected by the noise, this procedure ensures equal weigh- 
tage for presence of noise at all scales. Reconstructing the time series signal 
with this perturbed set of wavelet coefficients gave us the noisy signal to be 
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cleansed. Our studies in this fashion did show that noise and signal separa- 
tion was achieved. For the subsequent studies we followed the usual way of 
corrupting the signal by additive noise, i.e., 

x(U) =x°(t i )+ri(t i ) , (6) 

where r]{tj) G [—.5, .5] is the noise with zero mean and uniform distribution, 
and J the number of available dyadic scales. We have taken data size of 16384 
(= 2 14 ) points for all these three systems, and so J = 14. 

In Fig. 1 we plot our observations for the Lorenz system. Fig. 1 (a) shows 
the power at different scales in the pure signal x°(ti). In Fig. 1 (b) we plot the 
scalewise power distribution after numerically differentiating the pure signal. 
We see that almost the entire power of the differentiated data is accumulated 
within the dyadic scales 4 and 9 (the signal power between scales 10 and 14 has 
disappeared by the process of differentiation). Fig. 1 (c) plots the scalewise 
power in the noisy signal x{U) when the pure signal is infected with noise 
r](ti) of a typical strength of 5% of the signal (that is, 5% of the difference in 
the maximum and minimum values of x°{ti)). Because of the relative larger 
contribution at all scales from the pure signal, compared to that from noise, 
it is impossible to distinguish between the two components, and the figure 
looks qualitatively very similar to Fig. 1 (a). However, when we plot the 
scalewise power distribution of the differentiated noisy data in Fig. 1 (d), the 
signal contribution can easily be identified and also compared with the plot 
in Fig. 1 (b). It is to be noted that the difference in the values of the two 
peaks in Figs. 1 (b) and (d) arises because of the power being normalized by 
the respective total signal power. The contribution due to noise shows up in 
the finer scales. A clear minimum with close to zero value separates out two 
distinct regions. Fig. 1 (e) exhibits a small segment of the signal after the 
noise has been successfully removed following the procedure outlined above. 
All the three signals - pure, noisy, and cleansed - are overlaid for the sake of 
comparison. 

In Fig. 2 we show the results for the Autocatalator and the Rossler reacting 
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systems. Fig. 2 (a) and (c) plot, respectively, the scalewise power distribution 
for noise-infected signals obtained from the two systems. Like in the case of 
Lorenz system, it is evident that here also one cannot distinguish the noise 
and signal components. Fig. 2 (b) and (d) exhibit scalewise power profile for 
the differentiated data of the two signals respectively. The clear separation is 
again obvious. 

In order to quantitatively estimate the efficiency of our denoising method, 
we have made the following error estimation (Kostelich and Schreiber, 1993) 
for the above three model systems. Since in all these cases the pure signal 
is known, a measure of the amount of error present in the cleaned data is 
obtained by taking rms deviation of the cleaned signal x(ti) from the pure 
signal x°(ti) as follows, 

1 N 1 ,„ 

E={jr y E(Ht i )-x°(t i )) 2 ) , (7) 
iv 1=1 

where N is the length of the time series. Similar quantity E for the noisy 
data x(ti) is also computed. The condition E/E < 1 guarantees that noise has 
been successfully reduced. The error estimator E/E is a natural measure for 
noise reduction when the original dynamics is known (Kostelich and Schreiber, 
1993). In Fig. 3 we plot E/E against noise strength, for the three model 
systems. We see that for the entire range of noise values, and even with noise 
level as high as 10% of the signal exhibiting chaotic dynamics, E/E remains 
appreciably below unity. Thus the plot demonstrates the efficiency of the 
approach. Different wavelet basis functions may change the nature and also 
improve the efficiency further. 

We now discuss our method when applied to raw data obtained from two 
real chemical systems. In the first system, the time series data was obtained 
from the measurements of the pressure fluctuations in a fluidized bed, which 
consists of a vertical chamber inside of which a bed of solid particles is sup- 
ported by an upwardly moving gas. Our system used a bed of silica sand 
particles (of mean diameter 200 microns) with a settled height of 500 mm, 
fluidized by ambient air in a transparent vessel 430 mm across and 15 mm 
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wide. Beyond a critical inlet gas velocity, viz. the minimum bubbling velocity, 
the gas passes through the bed in the form of bubbles, thereby churning the 
solid and gas mixture in a turbulent manner. The time series data have been 
taken by measuring the pressure fluctuations inside this mixture, relative to 
atmospheric pressure, using a pressure transducer attached to a probe inserted 
into the fluid bed. The bed was operated at an inlet gas velocity of 0.85 m/sec, 
and the pressure fluctuations were recorded at a sampling rate of 333 Hz (333 
data points per second). As a standard procedure, we normalize the data by 
subtracting mean and dividing by standard deviation (Constantinides, 1987; 
Bai et. al., 1997). In Fig. 4 we show the results obtained after the data have 
been subjected to denoising. Fig. 4 (a) shows the power distribution at differ- 
ent scales in the original experimental signal, while in Fig. 4 (b) we plot the 
scalewise power profile of the differentiated data. Again one clearly sees the 
two distinct contributions due to the noise and signal components. Fig. 4 (c) 
shows short segments of the denoised signal and the original signal which is 
overlaid for comparison. The cleaned signal is seen to be smooth indicating 
that the noise has been removed. 

In the second chemical system, the time series data was obtained by sam- 
pling a measure related to the conductivity in a 3 liter liquid surfactant man- 
ufacturing experiment, at a sampling rate of 500 Hz. The time series is highly 
nonstationary since at various stages the operational parameters are altered 
(increasing the temperature for certain duration, then adding actives to the 
liquid, etc.). We studied unfiltered noisy data sets from the experiments, in 
order to check if our method can filter the noise out and also bring forth some 
intrinsic features of the system. We used our denoising algorithm to treat this 
data set in a slightly different way. The aim was to remove the finer scales 
from the differentiated data one by one, starting from the lowest (dyadic scale 
1) and gradually going up, so that at each stage (after integrating the data) 
the observable frequencies in the filtered signal may be related to identifiable 
physical sources. Fig. 5 (a) shows a small segment (1 second long) of the 
noisy data. In Fig. 5 (b) we plot, on the same scale as in the earlier figure, 
the filtered data, using our method to remove the lowest dyadic scale 1. One 
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can now clearly identify a 50 Hz component, due to the signal from electrical 
power supply (the 'net frequency'). By removing scale 2 alongwith scale 1, 
the net frequency goes away, and the filtered data exhibits a 13 Hz component 
superimposed with occasional spikes. This 13 Hz signal shows up clearly in 
the filtered data with scale 3 also removed. The same Fig. 5 (b) shows this 
data, overlaid on the data with scale 1 removed. This 13 Hz may have arisen 
from the stirring device which has two blades and revolves with 260 rpm, core- 
sponding to approximately 10 Hz. The electronic signal had an antialiasing 
feature of no more than 250 Hz and therefore aliasing (beating) may be ruled 
out. It may also be mentioned here that the Fourier power spectrum of the 
denoised signal shows a spike at 50 Hz frequency, whose removal resulted in a 
residual spectrum consisting mainly of a background continuum without any 
appreciable peak around 13 Hz. This study with the present example sug- 
gests that the wavelet transform methodology offers considerable benefits in 
the recovery of intrinsic signal components. 

Summary 

We have presented a new and alternative algorithm for noise reduction using 
discrete wavelet transform. We believe that our algorithm will be beneficial 
in various noise reduction applications, and that it shows promise in develop- 
ing techniques which can resolve an observed signal into its various intrinsic 
components. In our method the threshold for reducing noise comes out au- 
tomatically. The algorithm has been applied to three model flow systems - 
Lorenz, Autocatalator, and Rossler systems - all evolving chaotically. The 
method is seen to work quite well for a wide range of noise strengths, even 
as large as 10% of the signal level. We have also applied the method suc- 
cessfully to noisy time series data obtained from the measurement of pressure 
fluctuations in a fluidized bed, and also to that obtained by conductivity mea- 
surement in a liquid surfactant experiment. In all the illustrations we have 
been able to observe that there is a clean separation in the frequencies covered 
by the differentiated signal and white noise. However, if the noise is colored, 
a certain degree of overlap between the signal and noise may exist even after 



10 



differentiation. For this complex situation, the method needs to be improved 
upon. 
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Figure Captions 



Fig. 1. Plots for Lorenz system, with parameter values as stated in Table |l|. 
16384 (= 2 14 ) data points are considered. Scalewise power distribution is 
plotted against the dyadic scale, (a) for the pure signal (the interpolated 
line through the data points is drawn for visualization), (b) using the 
data obtained after differentiating the signal, (c) for the signal corrupted 
by noise, and (d) using differentiated data of noisy signal. A segment 
of cleansed signal is shown in (e), alongwith the pure and noisy signals 
overlaid for comparison. 

Fig. 2. Plots for Autocatalator and Rossler systems, in (a), (b) and (c), (d) re- 
spectively, for parameter values as in Table [l]. Scalewise power profile 
plotted, (a) for the noisy autocatalytic signal, (b) using data after dif- 
ferentiating the signal, (c) for noisy Rossler signal, and (d) using the 
differentiated data of noisy signal. 

Fig. 3. The error estimator E/E plotted against the noise strength for all the 
three systems. 

Fig. 4. Plots for the fluidized bed experiment, (a) Scalewise power profile is 
shown, for (a) the experimental signal, and (b) the data after the signal 
has been numerically differentiated. A small segment of the cleansed 
signal is shown in (c), alongwith the original signal for comparison. 

Fig. 5. Plots for the liquid surfactant experiment, (a) A segment of the original 
noisy data, 1 second long, (b) The filtered data, with scale 1 removed, 
and with scales 1, 2 and 3 removed (on the same axes-scales as (a)). 
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Table 1: Various model systems studied, alongwith their parameter values and 
nature of dynamics. 



Systems 


Lorenz 


Autocatalator 


Rossler 


Dynamical 
equations 


dx/dt = —ax + ay, 
dy/dt = Rx — y — xz, 
dz/dt = —bz + xy. 


dx/dt = 1 — x — Da\xz 2 , 
dy/dt = (3 — y — Da 2 yz 2 , 
dz/dt = 1 - (1 + Da 3 )z 
+a(Daix + Da 2 y)z 2 . 


dx/dt = —y — z, 
dy/dt = x + ay, 
dz/dt = b 

+z(x — c). 


Parameters 
chosen 


a = 10, R = 28, 
b = 8/3. 


D ai = 18000, Da 2 = 400, 
Da 3 = 80, (3 = 2.93, 
a = 1.5. 


a = .398, 6 = 2, 
c = 4. 


Dynamics 


Chaotic 


Chaotic 


Chaotic 
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